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Q . ABSTRACT 



Recent observations of quasi-periodic oscillations in the aftermath of giant flares in 
soft gamma-ray repeaters suggest a close coupling between the seismic motion of 
the crust after a major quake and the modes of oscillations in a magnetar. In this 
£N| ■ paper we consider the purely elastic modes of oscillation in the crust of a neutron 

star in the relativistic Cowling approximation (disregarding any magnetic field). We 
determine the axial crust modes for a large set of stellar models, using a state-of-the- 
art crust equation of state and a wide range of core masses and radii. We also devise 
■ useful approximate formulae for the mode-frequencies. We show that the relative crust 

thickness is well described by a function of the compactness of the star and a parameter 
describing the compressibility of the crust only. Considering the observational data 
for SGR 1900+14 and SGR 1806—20, we demonstrate how our results can be used to 
. constrain the mass and radius of an oscillating neutron star. 
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1 INTRODUCTION 

Asteroseismology aims to probe stellar physics via various observed modes of vibration. To be successful in this endeavour one 
» I ■ needs both accurate observations and detailed theoretical models to test the observations against. In the case of neutron stars 
we have until very recently had neither reliable observations nor detailed theoretical models. The situation appears to have 
changed with the observations of quasiperiodic oscillations (QPOs) following giant flar es in three soft gamma-ray repeaters 
(SGRs) ijlsrael et, alJl200.4 Istrohmaver fc Wattell200,4 EuOrl IWatts fc St,rohmaverll200fif) . Analysis of the relevant X-ray data 
has unveiled a number of periodicities , with frequencies that agree reasonably well with the expected torsional modes of the 
neutron star crust [e.g. JPuncarj |l998'll . These observations are tremendously exciting because they will allow us, for the very 
first time, to test our neutron star oscillation models. They also provide strong motivation to improve these models. After all, 
in order to be accurate enough to make a comparison to the observations meaningful, a model must be fully relativistic. It 
should also allow for the presence of a strong magnetic field and possible superfluid components. 

In this paper we take some modest steps towards more realistic models of pulsating neutron stars by including the crust 
elasticity in a relativistic calculation of axial oscillations. To simplify the problem we work within the Cowling approximation, 
i.e. we neglect the dynamical nature of spacetime. This is expected to be an accurate approximation for axial modes in the 
crust since they originate in a region of relatively low density and should not lead to sig nificant variation s in the gravitational 
field. We are not considering the class of axial gravitational- wave w- modes at this point llKokkotaslll994l) . apart from formally 
noting how they would appear in our formulation (see Appendix A). We calculate axial modes for an up-to-date crust equation 
of state (not very controversial) for a large set of core masses and radii, corresponding to different (unspecified) supranuclear 
equations of state in the core and (similarly unspecified) central pressures. Our approach allows us to largely ignore the 
detailed structure of the core which is only poorly constrained by observations and theoretical considerations. We show how, 
when used in conjunction with observed mode frequencies, our numerical results can be used to put constraints on the global 
equation of state. We also determine an approximate solution to the mode-problem. This leads to useful formulae that show 
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how the axial modes depend on the key parameters, the star's mass and radius, the crust compressibility and the shear 
modulus. When compared to our numerical data these approximations are surprisingly good. We show how they can be used 
to deduce the main stellar parameters from a set of observed frequencies, thus outlining a complete asteroseismology analysis. 

Our study does not at this point account for the presence of a magnetic field. This obviously means that one should 
be careful before assuming that our results can be used to interpret the magnetar observations. However, we feel that the 
magnetic problem is still somewhat beyond reach. The main reason for this is the strong cou pling bet ween any mot ion in 
the crust and Alfven waves in the core. As we have recently argued iciampedakis et alJ 20061) [see also iLevml teOQctl l. this 
coupling is likely to lead to global oscillations which involve significant core motion. This problem is computationally orders 
of magnitude more difficult because the magnetic field plays an active role. Having said that, it may well turn out to be the 
case that the modes that dominate the observed signal remain quite close to the pure crust modes, desp ite the presenc e of the 
magnetic field. One can argue that this should be the ca se if the modes are excited by crust "cracking" founcanll998l) . where 
the bulk of the energy is deposited into elastic motion Jciampedakis et alj l2006). This is, of course, a qualitative argument 
and better calculations are needed to establish whether or not it is quantitatively useful. 

Before we proceed, it is worth noting that the approach to neutron star asteroseismology described in this paper is directly 
applicable to weaker magnetic field pulsars provided that they are not too rapidly rotating. Glitches in these objects may 
excite torsional oscillations in the crust, even though such oscillations are yet to be detected. 



2 PERTURBATION EQUATIONS 

Our study is motivated by a key question for both theorists and observers: To what extent can we use the observed QPOs 
to probe neutron star physics? We attempt to answer (at least partially) this question by calculating axial oscillation 
modes for a realistic neutron star crust using a fully relativist i c theory of elasticity that was deve loped in a recent series 
of pape rs iKarlovini fc Samuelssonl 1200.4 iKarlovini et alJ 120041 : iKarlovini fc Samuelsson! I2004L l200fih based on the work of 
ICarter fc Quintanal (jl972T ). Since we want to emphasise the results and the possible implications for asteroseismology rather 
than various computational technicalities we will consider a simple neutron star model that (we believe) captures most of the 
key features that govern torsional crust oscillations. 

In our model the core has radius R c and mass M c but we keep the equation of state (EOS) unspecified. We can do this 
since, wi thin our approxima tions, there is no coupling between the crust and the core fluid. The EOS o f the crust is described in 
detail bv lSamuelssonl l)2003|) and its essen tial parts are based on the work of lHaensel fc Pichonl l|l994F) and lDouchin fc Haensell 
The shear modulus is taken from Oeata & Ichimaru (1990) and we assume that the temperature is zero. 

We let the background star be static and spherically symmetric. This means that the spacetime metric is given by 

ds 2 = ~e 2v At 2 + e 2X dr 2 + r 2 (d6 2 + sin 2 <9d^ 2 ). (1) 

The perturbations are treated in the Cowling approximation, i.e. we neglect metric perturbations. This should be a good 
approximation provided that the oscillation modes are confined to the relatively low density region represented by the crust. 
The perturbation equations for the crust, assuming perfectly elastic matter, are derived in Appendix^ These equations allow 
for the solid to be anisotropic. The anisotropy manifests itself by the existence of two different speeds of shear waves, denoted 
by v r (radially propagating with polarisation in an angular direction) and vt (propagation in an angular direction, polarised 
in the mutually orthogonal angular direction). Throughout the paper we will use a subscript r to mean radial and t to mean 
tangential (to the spherical symmetry surfaces). Similarly, there will be two different pressures, p r and p t . In the isotropic 
limit we denote the quantities without subscripts, i.e. v = v r — vt and p = p r — pt- 

As a further simplification, we assume that the core is unable to support traction. This effectively isolates the crust, since 
we neglect, e.g. viscosity and magnetic fields. While the for mer should be a reasonab le approximation, the latter is not valid 
in the context of magnetars. As we have argued elsewhere jGlampedakis et alJl200(3) . the crust-core coupling time scale in a 
magnetar is comparable to the oscillation period of the crustal modes which implies that one ought to consider global MHD 
modes. As this proble m is seriously challenging, we prefer to focus on the non-magnetised case here. After all, if we take the 
toy-model discussed bv lGlampedakis et alJ <|2006F) at face value then it seems plausible that the modes that are mainly excited 
by crust cracking in the magnetic problem may be close to the pure crust modes determined in the non-magnetic case. 

In Appendix [X] we show that the axial perturbation equations for an elastic solid in the Cowling approximation can be 
written 

F" + A'F' + BF = 0, (2) 
where a prime denotes differentiation with respect to the Schwarzschild radius r and 

e A = r 4 e v ~ x (p + p t )vl, (3) 
v 2 (l -!)(/ + 2)" ([) 



B = 
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Here F describes the amplitude of the fluid oscillations, p is the energy density, the integer I is the usual angular separation 
constant that enters when the perturbation variables are expanded in spherical harmonics Y l m (8,4>) and u> is the angular 
frequency. 

In the case that we are considering, the boundary conditions require that the traction vanishes at the top and the bottom 
of the crust. In terms of our independent variable F this leads to 

T = e A F' = at r = R c and r = R. (5) 

It is convenient to use T to reduce to a first order system, 

T' = -e A BF, (6) 

F' = e' A T. (7) 

This system is suitable for numerical integration since there are no derivatives of the equation of state parameters, e.g. the 
shear modulus or the density ( which are only known i n tabular form). In the isotropic limit these perturbation equations are 
equivalent to e.g. those used bv lYoshida fc Leel l|2002F) . 

We have used these equations to calculate axial crust modes for a large number of stellar models with different core 
masses and radii, implicitly corresponding to a variety of supranuclear equations of state in the neutron star core. We want to 
use these results to develop a workable strategy for asteroseismology. That is, we would like to i) identify the key parameters 
that govern the various modes, and ii) try to represent the results in such a way that a parameter "inversion" becomes 
possible given actual observations. The basic idea behind this kin d of analysis has already bee n discussed for neutron star 
oscillations that radiate gravitational waves at a significant level llAndersson fc Kokkotasl 19981) . Since the analysis is much 



aided by "analytic" formulae, we will develop a suitably simple approximation scheme for the axial modes. This will lead to 
a set of empirical expressions that can be used to address the inverse problem. 



3 SIMPLE ANALYTIC APPROXIMATIONS 

We want to find an approximation to the axial crust modes. To do this, we consider Eq. and introduce a new Regge- Wheeler 
type radial coordinate x through 

^=e- A (8) 
dr 

The perturbation equation then becomes 

^+e 2A BF = 0, (9) 
dx 2 

i.e., it is written on a form that lends itself to a WKB-type approximation. Since we are interested in making progress 
analytically we consider only the lowest order approximation. This means that we assume that the solution can be written 

F = Cxe imx) +C 2 e-* W{X \ W(x) = f e A B 1/2 dx = f B 1/2 dr. (10) 

J R a J R a 

At the base of the crust (r = R c ) we need to ensure the vanishing of the traction. Hence, we impose the boundary condition 

F' = iB 1/2 (d-C 2 ) =0 => Ci=C 2 . (11) 
Meanwhile, the analogous condition at the surface (r = R) implies that 

F' = iCiB 1/2 (e iW(fl) -e- w(fl) )=0 => W(R) = [" B 1/2 dr = utv. (12) 

To make further progress, we next make the approximation that the shear speeds are constant (an approximation which 
is good throughout much of the crust) and that v and A are constant (which is a reasonable approximation since the crust 
mass in negligible compared to that of the core). Then, assuming that 

a a „ 2 + 2) 
u >e vt 5 , (13) 



we may Taylor expand _B 1//2 to get 

(14) 



,1/2 X-v Id 
V r 



e 2v v 2 (l-l)(l + 2) 
2u 2 r 2 



which may be integrated to yield, using Eq. 1121 . 



A 2RR C ' v ' 

where A = R — R c . This equation provides a useful first approximation to the frequencies of the axial crust modes. 
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Figure 1. In the left panel we show that the ratio of the frequencies of the fundamental modes for I = 2 and I = 3 agrees well with 
that predicted by our (I — l)(l + 2) scaling which in this case gives a ratio of \/5/2 0.6324555 (shown as a solid line). The small 
deviation seen for low /3 = M/R is typical for the first few Vs. The remarkable agreement for typical neutron star values j3 ~ 0.2 should 
be noted. It is also worth remarking that our numerical code has an accuracy of a few times 10 — 5 , which means that some of the 
remaining spread in this range may be due to numerical error. In the right panel we show that the numerically computed ratio of the 
frequencies of the fundamental modes agrees well with that predicted by our (I — l)(l + 2) scaling (shown on the left as solid lines). 
The 1(1 + 1) prediction is show on the right as dashed lines. Of course, for large I (so that the plane wave approximation is valid) we 
get (I - l)(l + 2) = I 2 + I - 2 -> 1(1 + 1). Nevertheless, the di screpancy for the first few values of £ may lead to an erroneous multipolc 
assignation for a given observed mode-signal. 



Before considering the overtones, let us discuss the case of the fundamental crust modes, which correspond to n = 0. For 
this case we immediately find 

^ ^ { i-m + 2) (n=o) (l6) 

However, this result shows that our approximation is not consistent for the fundamental mode. It is clear that the condition 
1131 1 is not satisfied throughout the crust. Nevertheless, our approximate resu lt provides useful insights into the scaling with 
various parameters. We can compare it to, for example, the formula used bv IPirol ^2005^1 [his Eq. (9)]. This is the standard 
estimate, which is arrived at via a plane-wave approximation, and it can be written 

^^|^, (17) 

where v is the isotropic speed of shear waves. This result differs from ours in two important ways. First of all, the I dependence 
is different. Secondly, our formula replaces R 2 by RR C . 

It is easy to demonstrate that the difference in the /-dependence is important, and that our formula represents the correct 
behaviour. In order to do this, let us consider the ratio between the fundamental quadrupole (/ = 2) mode frequency and 
various higher multipoles. In the ratio the dependence on, for example, the shear wave speed disappears and we are left with 
a simple scaling with I. If we also compare the approximate formula to the results obtained from solving the full problem 
numerically we get insight into the accuracy of some of our assumptions, in particular whether it is reasonable to treat the 
shear wave speed as if it were constant and effectively "the same" for the different modes. This comparison is illustrated in 
figure The left panel shows the ratio between the fundamental I — 2 and I — 3 modes for our approximate formula and 
our total sample of 9000 numerical models with parameters in the range 8 < R c < 16 km and 0.8 < M c /Mq < 2.3. From 
the figure it is clear that the numerical results scale with I in the way suggested by our (inconsistent) approximation. This is 
further demonstrated by the data in the right panel of figureQ which extends the comparison to higher multipoles. This figure 
shows that an assumed scaling with 1(1 + 1) is likely to lead to an erroneous identification of the various higher multipoles 
that may be present in an observed signal. It may, fo r instance, be noted in table U that we ar rive at a different multipole 
identification compared to Istrohmaver fc Watts] (120051) for SGR 1900+14 and llsraelet all J2005I) for SGR 1806-20. 

Let us now consider the various overtones n ^ for any given I. We then need to solve the quadratic equation 115H . 
Expanding the resulting square root in the small parameter A/R and ignoring the negative root we get 

u-xn-KVr (. , 2a (Z~1)(Z + 2K 2 A 2 1\ 
" 6 —{ 1 + e 2^ ^RR-^J- (18) 

The second term in the parenthesis is negligible for moderate Z. It is worth noting that it is v t that affects the n — modes 
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whereas the n > modes are primarily determined by v r . It should also be emphasised that condition 113H holds for the 
overtones, which means that our approximation is consistent in this case. 

In order to estimate the overtone frequencies for any given M and R (say) we need to provide the crust thickness A(R, M) 
which then leads to R c . In Appendix [5] we show that the crust thickness is well approximated by 

i- (§«"-)"'• M 

where /3 = M/R is the stellar compactness and a is a parameter that depends on the equation of state and which essentially 
measures some average compressibility of the crust. As discussed in Appendix B, the relevant value for our crust equation of 
state is a = 0.02326. 

We want to combine these various approximations to get explicit expressions for the mode frequencies that can be used to 
compare to either numerical results for detailed relativistic models or observational data. To do this, we first assume that the 
crust is isotropic by setting v r = v t = v. If we then compare Eq. JTU and JTHJ to our full numerical results we find that they 
provide reasonable approximations for the same value of the shear wave speed v, provided that we adjust the factor of 1/2 in 
1161 . Since we already knew that Eq. 1161 was not derived in a consistent way, this manipulation does not seem unreasonable. 

The shear speed varies very little in the crust and has a value ~ 10 s cm/s. Writing v = vq x 10 s cm/s, where vo is a 
dimensionless parameter we may write our frequency estimates as 



/~11S5 Vo V0- 1 )0 + 2 ) / (1-705 - 0.705/3.) (0.1055 + 0.8945/3.) 
-Rio 2 y P* 

/^ 473.77n-^_(0.1055 + 0.8945/3.), (n > 0) (21) 
-Rio 

where .Rio = R/10 km and /3* = /3/0.2068 (/3. = 1 for a star with radius 10 km and mass 1.4Mq, obtained using the most 
recent value for G etc.). Comparing to our numerical data we find that Vo ~ 2.34. This value may appear to be surprisingly 
large given that, for our EOS, the maximum shear speed is about 2 x 10 8 cm/s. However, it should be remembered that 
the parameter vo is obtained through comparison with the full numerical solutions and therefore represent some "weighted 
average" shear speed. In this (implicit) averaging geometrical factors and explicit EOS dependence may well contribute factors 
of order unity. In a sense, the replacement of R with the geometrical mean radius of the crust V RRc represent the first order 
geometrical effect. It is also worth pointing out that the simple linear scaling with n does not represent the behaviour found 
in the numerical study. Again, this is not surprising as the eigenf unctions will have different forms and hence depend more 
strongly on the EOS in different parts of the crust. In general the mode frequencies seem to grow less rapidly than linear as a 
function of n. Another interesting point to note is the absence of an overall redshift factor in the estimated overtone frequencies. 
This somewhat counter-intuitive result arises since the redshift factor in l|18|l is essentially cancelled by the corresponding 
factor in the formula for the crust thickness. 

The accuracy of our approximate formulae is illustrated in Figure|5| When comparing the approximations to the numerical 
data we note that, for the fundamental modes, the error is a smooth function of 0. The situation therefore lends itself to 
improvement by fitting. Playing the same game with the overtones we find 



f~ 27 65 1 VQ-l)(Z + 2) / (1.705 - 0.705/3. ) (0.1055 + 0.8945/3.) ft 

1 ~ ' R 10 2 y /3. 1.0331/3. - 0.0331 ' 1 ' 1 ' 

f w 1107.3^(0.1055 + 0.8945/3.) — , (n > 0) (23) 

These frequency estimates are evaluated against numerical data in figure [5] We see that they perform extremely well, with 
relative errors typically smaller than 1%. Clearly, our various approximate formulae provide a very accurate representation of 
the full numerical results. We therefore anticipate that they will prove useful for any attempts to infer physical parameters 
from observational data. 



4 SEISMOLOGY 

In this section we discuss how one can use our numerical and approximate results for axial modes in the neutron star crust 
to infer physical parameters from observations, eg. the QPOs found in the tails of magnetar flares. To carry out this analysis 
we assume that the observed QPOs in the frequency range ~ 28 — 160 Hz are connected with fundamental (n = 0) elastic 
modes with different l's, whereas the QPO at ~ 626 Hz is an overtone (see below for a brief summary of the observations). 
It is important to note that these assumptions may not be valid, especially since one can argue t hat the crust motion excites 
motion also in the core fluid via the magnetic field coupling ijciampedakis et al.l2006llLevinl2006l) . If this is the case, then the 
mode problem that we have solved in this paper does not provide a complete picture. However, even though our analysis may 
not be entirely realistic in this sense we believe that our discussion provides a useful example of the strategy. We should also 
stress that, until fully relativistic modes for magnetised stars are calculated, this analysis may be the best one can expect. 
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Figure 2. In the left panel we compare various estimates for the fundamental I = 2 torsional mode. We show the relative error, compared 
to numerical data, of the frequency estimates as functions of /3. In each case we have rescaled the estimates to give approximately the 
numerically computed value of the frequency at /3* = 1 (see the main text). This corresponds to choosing an appropriate value of the 
shear speed. Comparing Piro's result (c) with ours (a) we see that the (RR C )~ 1 ^ 2 scaling performs better than the standard R —1 . We also 
note that the extrapolated "general relativity correction factor" of Duncan (d) not only makes the estimate less accurate, but actually 
destroys the clean /3-dependence. Also shown is our fit of the frequencies (b) which is accurate to within 0.2 %. In the right panel we 
show the relative error in the estimate for the frequency of the first overtone both with (b) and without (a) the empirical fitting factor. 

4.1 Observations 

Before considering the inverse problem of finding global neutron star properties from observed spectra let us briefly discuss 
the observations. So far there have been three recorded giant flares, all associated with SGRs. 

SGR 0526-66: This flare took place in March 5 1979 jMazets et alJll97Sh . Only a weak QPO was detected at w 44.5 Hz 
jBarat et al.lll983fo . Due to the meagre data and the fact that this feature was detected in the peak of flare with significant 
dead-time effects in the data we will not consider this observation further. 

SGR 1806-20: Erupted December 27 2004 Jfiurlev et al.l2005l:IPalmerl2005l) . llsrael et alJ(l2005l) detected QP Os at 18.1±0.3, 
30.4±0.3 and 92.5 Hz in the data from the Rossi X-ray Timing Explorer fRXTE ). \Watts fe StrohmaveJ (2006) later examined 
the data from the Ramaty High Energy S olar Spectroscopic Imager (R HESSI) and found QPOs at 17.9 ± 0.1, 25.7 ± 0.1, 30 
(weaker), 92.7 ±0.1 and 626.46 ± 0.02 Hz. Istrohmaver fc Watt! I2OO6I) have recently detected a feature at » 150 Hz 1 as well 
as showing that the w 30 Hz QPO is better fitted by 29 Hz. In our analysis we exclude the lower frequency QPOs. This choice 
is guided by our recent toy problem ijoiampedakis et, alJl200f^l . Hence, we consider the 29, 92.7, 150.3 and 626 Hz modes. The 
inclusion of the 150 Hz QPO does not change the conclusions substantially (which, incidentally, support the identification of 
the 29 Hz QPO with the I = 2 mode and that our scaling with I is the correct one). In table Q we display these frequencies 
together with the modes that we identify them with. 

SGR 1900+14: A giant fl are from this magnetar was detected on August 27 1998 jHurlev et a,lJll99fll IFeroci et, alJll99fll ). 
Istrohmaver fc Watts! (|2005i) found evidence for QPOs at 28 ± 0.5, 53.5 ± 0.5, 84, 155.1 ± 0.2 Hz in the data from RXTE. In 
our analysis we use all of these frequencies. The modes that we identify them with are listed in tabled 

4.2 Seismology with the analytic expressions 

It is straightforward to use the analytic approximations that we derived in the previous section to analyse the observed 
frequencies and extract the key neutron star parameters. 



1 Thi s feature was brought to our intention via private communication with Watts. After this work was completed Strohmav er fc WattJ 
(2006) published their results from a reanalysis of RXTE data. In addition to the new 150 Hz feature several higher frequency QPOs 
were found. These have not been considered in this work. 
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SGR 1806- 


-20 


SGR 1900- 


1-14 


f [Hz] 


mode 


f [Hz] 


mode 


29 




28 ±0.5 


0*2 


92.7 ±0.1 


0*6 


53.5 ±0.5 


0*4 


150.3 


0*10 


84 


0*6 


626.46 ± 0.02 


1*1 


155.1 ±0.2 


0*11 



Table 1. Observed QPO frequencies and suggested corresponding elastic crust modes. The identified modes are denoted according to 
n ti. Note that the true ^-dependence differs from the standard 1(1 + 1) which explains why our identifications are different from those of 
Strohmayer & Watts (2005) (SGR 1900+14) and Israel et al. (2005) (SGR 1806-20). It may also be noted that the rather large allowed 
error (6 %) in our analysis makes it possible to identify the 150 Hz QPO in SGR 1806—20 with both an / = 9 and and / = 10 mode in 
different regions of the M-R parameter space. However, the presence of the overtone at 626 Hz marginally rules out the 1 = 9 option. 



Let us first assume that we observe the fundamental quadrupole mode together with the first overtone (n = 1). In the data 
for SGR 1806—20 it seems reasonable to assume that the first is represented by the 29 Hz QPO, while the latter corresponds 
to the 626 Hz mode. From our approximate formulae we immediately see that the ratio of these modes provides an expression 
that only contains the compactness /?*. In effect, the data provides a curve in the M — R plane on which the true stellar model 
should lie. Solving this constraint for the compactness we find that j3 ~ 0.12, i.e. R ~ 8.1M. From Eq. 1191 we also see that 
the relative crust thickness is A/R ~ 0.17. We can then insert this value for (3 in the expression for (say) the fundamental 
mode. This allows us to solve for the radius, and we find that R w 11.4 km, which means that M « 1.41 km « O.96M0 and 
A « 1.9 km. These values do not seem unreasonable, even though the inferred mass is quite low and the crust surprisingly 
thick. Most likely there are systematic effects due to the magnetic field. These are, of course, not accounted for in our analysis. 

Next consider the fundamental modes for different values of /. First note that the scaling with I allows us to immediately 
work out the ratio of the different mode frequencies. Once we assign the ~30 Hz QPO to the fundamental I = 2 mode we can 
infer the multipoles of the various higher frequency modes in the observed data. The conclusions of this analysis are given in 
Tabled Once we have determined the various multipoles, our approximate formula for the fundamental modes can be used 
to constrain the stellar parameters. Let us take the quadrupole mode as an example. If the frequency is known, then the 
approximate formula gives the stellar radius as a function of the compactness f3. Since Al = /3R we again have a constraint 
curve in the M — R plane. In order to be consistent with the results for the n = 1 overtone, the two curves must intersect. 
The point of intersection immediately provides us with the mass and radius of the star. Available results for higher i-modes 
provide further constraints that can be used to verify the consistency /accuracy of the inferred stellar parameters. Of course, 
since our model is an idealisation (a spherical star with an isotropic crust etcetera) one would expect the use of real data to 
lead to a spread of the various curves in the M — R plane. This uncertainty can to some extent be used to assess the reliability 
of the parameter extraction. 



4.3 Using numerical results to identify parameters 

In our numerical code we solve the background Einstein equations together with the perturbation equations, specifying the 
core mass and radius. In the absence of discontinuous tangential pressure at the crust-core interface neither this boundary 
nor the top of the crust present a ny difficulties. [If there i s a ju mp discontinuity in t he ta ngential pressure we have one 
more parameter in the problem, see lKarlovini fc Samuelssonl J2003 ) and lKarlovini et alJ J2004t) for discussions.] It is therefore 
straightforward to integrate the equations both inwards and outwards to a matching point where we make sure that the 
Wronskian of the system (|SJ and Q vanishes. We have checked that the results do not depend on the choice of matching 
pressure. The simple nature of our model allows us to quickly compute many modes (2 < I < 14, < n < 3) for a large set of 
core masses (45) and radii (200) numerically. 

Our aim is to use this data set to determine models that match all the observed frequencies listed in table (to some 
accuracy). As indicated above, one has to allow for some level of uncertainty when comparing the observed data to our 
numerical calculations. In our analysis we have, rather arbitrarily, chosen to allow a relative error of 6 % in the frequencies. 
For each of the two flares we associate a subset of the observed QPO frequencies with axial oscillation modes and search our 
computed data set for stellar models which has all the required frequencies, for some sequence of Vs. For SGR 1900±14 we 
use all four observed QPOs and infer that the excited modes correspond to I = 2, 4, 6, 11. For SG R 1806—20 we assume that 
the frequencies lower than about 30 Hz are not directly associated with axial crust oscillations [see loiampedakis et al.l {2006) 
for a possible explanation of these modes]. We are then left with a sequence of three fundamental modes which we identify 
with the multipoles I = 2,6, 10. The last one, when only compared to the other fundamental modes, could also be allowed to 
have I = 9 in parts of the parameter space. However, for SGR 1806—20 we also a higher frequency QPO which we associate 
with an n = 1 overtone with arbitrary (but low) I. This marginally restricts the higher fundamental mode to / = 10. 



8 Samuelsson & Andersson 



2.5 



M 

'o 

00 



23 1.5 




i 1 l 

10 12 14 

Radius [km] 



i 1 i 1 
SGR 1900+14 



s 



G300 



G240 



\ 



12 

Radius [km] 



14 



16 



Figure 3. Seismology analysis based on our numerical axial-mode results. For the flare in SGR 1806—20 (left panel) we associate the 
lower frequency QPOs (approximately 29, 93 & 150 Hz) with fundamental (n = 0) modes with I = 2, 6, 9 or 10 respectively. The inclusion 
of the 150 Hz mode does not change the picture substantially. The higher frequency QPO (626 Hz) is associated with a n = 1 mode of 
arbitrary I. The models allowed by the lower QPOs form a rather broad region that is orthogonal to the line corresponding to models 
that have an n = 1 overtone of the right magnitude. As discussed in the main text, the true stellar parameters should belong to the 
region of overlap. For the flare in SGR 1900+14 (right panel) we use all the reported QPOs (approximately 28, 53, 84 & 155 Hz) and find 
that the sequence of Z's is (2,4,6,11). The larger number of modes leads to a narrower allowed region for the given tolerance. However, 
the lack of overtones in the data means that the stellar parameters are much less constrained in this case. 



The results of the analysis are displayed in Figure In each panel we plot the mass-radius relation (as li nes) for a few 
propo sed EOSs. The equations of state are chosen for illustration only and include the A18-fo+UIX* mode l of lAkmal et alJ 
lll99Sl) with and with out a deconfined quark core (APR fc APRQ), two examples from iGlendennin j feOQCl) (G240 & G300) 
and the FPS EOS of IPandharipande fc Ravenhalll (^L989) (FPS). On top of this we display, as grey diamonds, the allowed 
masses and radii if the given set of observed low frequency QPOs are associated with elastic oscillations. In the plot of SGR 
1806—20 the models which have an n — 1 overtone is shown as dark squares. 

The numerical results nicely illustrate the strategy for neutron star seismology discussed for the approximate formulae. 
In particular, the case of SGR 1806—20 indicates that the detection of overtone modes is extremely valuable. Given both the 
fundamental mode and an overtone the stellar parameters can be constrained to a relatively narrow region of the M-R plane. 
This, of course, then provides useful constraints on the supranuclear core EOS. 



5 DISCUSSION 

In this paper, we have considered torsional oscillation modes in neutron star crusts. Using a state-of-art EOS in the crust 
but otherwise simplifying the problem as far as possible we have shown how global properties of the star may be deduced 
from observations of such modes. We have done this using both approximate formulae and a large sample of fully numerical 
mode-results. As an example we applied our methods to SGRs assuming that the post-flare QPOs are connected with pure 
crustal oscillations. The results show that the devised strategy for asteroseismology is, in principle, quite powerful. Of course, 
the basic model needs to be improved before the inferred stellar parameters should be taken too seriously. In particular, we 
expect systematic effects due to our neglect of the magnetic field, which obviously should be important for magnetars. One 
can immediately think of several ways to improve upon our calculations. Below we list the most important. 

For a magnetic star, one would not expect the QPOs to be exactly at the pure crustal frequencies. The toy model of 
iGlampedakis et alJ J2OO3) suggests that the QPOs are really global oscillations of the entire star, but that the modes that 
are predominately excited lie close to the pure crust modes. In addition one would expect the geometry and strength of the 
magnetic field to have a key role in determining which modes are excited and how strongly. Clearly, this question can only be 
settled in a realistic spheroidal model. Work on that problem should be strongly encouraged. 

The magnetic field will not simply lead to an offset of the mode frequency with respect to that o f a pu re crust mode. 
The nature of the modes will also be affected by the presence of the field. The considerations of IPirol ((2005) show that (in 
his plane-parallel geometry) the fundamental modes are only weakly affected by the magnetic field, whereas the overtones are 
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strongly dependent on the field. This is to be expected from our expressions 1221 and 123H where we see that the fundamental 
modes are determined by the tangential shear velocity vt while the overtones instead depend on the radial velocity v r . A 
magnetic field will boost the elasto-MHD speeds along the field direction and it therefore seems possible that the toroidal part 
of the field will affect the fundamental modes whereas the poloidal part will influence the overtones. Should this prove to be 
true for the spherical case it would be very interesting since it would offer a possibility of extracting information concerning 
the field geometry. That the geometry o f the magnetic fi eld configurat i on is essential to the behaviour of the torsional modes 
was emphasised bv lMessios et al] J200jl and recently bv lSotani et alJ teOQcj) . 

In our model we h ave assumed zero tempe rature. A non-zero temperature will lead to a systematic effect since the shear 
modulus will decrease JOeata fc Ichimarull 990). It will also mean that the crust will be slightly thinner as the ocean is deeper. 
In our calculations we have integrated the equations all the way to the density of iron at zero pressure. 

In the deep crust matter is not expected to be microscopically isotropic, but instead form so-called pasta phases. If these 
anisotropies extend to macroscopic (fluid-dynamic) scales, i.e. if they do not average out over larg e distances compared to th e 
lattice spacings, the radial and tangential shear speeds will be different, perhaps by a factor of two JPethick fc Potekhinll99gl) . 
Additionally, the crust may not be in a relaxed state, causing the background strain to influence the torsional oscillations. 
We have not taken this into account in our model. 

In the inner crust we expect that at least a portion of the dripped neutrons are in a superfluid state. It is, in fact, likely 
that a significant portion of the mass does not directly participate in the oscillations. Of course, the superfluid will n evertheless 
be set in motion via the effect of entrainment [for discussion and relevant references see lAndersson fc Comerl l|200fiF) ]. It would 
b e very interesting the examin e the effects of this mechanism on the mode spectrum using the recently developed formalism 
of lCarter fc Samuelssonl l)200fit) (which incidentally also allows for a MHD-type magnetic field) . 

Our models are static, i. e. non-rotating. Although this is a very good approximation for the slowly rotating magnetars it 
would be useful to extend our treatment if toroidal oscillations were to be observed in faster spinning pulsars. 

Our model is non-dissipative. In reality a number of processes, including the emission of gravitational and electromagnetic 
radiation, will damp the oscillations. Mutual friction between superfluid and normal components in the inner crust and viscosity 
will also be important. Observations of the QPOs suggest a damping time of a few tens of seconds in the SGRs. It would be 
interesting to estimate the characteristic damping times for each damping mechanism to see if they are consistent with the 
observations. In addition, if the emission mechanism of the observed X-ray signal was known, a lower limit of the amplitude 
of the motion of the crust could be derived. This might lead to an estimate of the maximum strain in the crust and hence an 
assessment of the viability of elastic motion as the origin of the QPOs. Clearly, if the maximum strain exceeds the breaking 
strain, the oscillations cannot be of elastic origin. 
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APPENDIX A: AXIAL PERTURBATIONS IN THE COWLING APPROXIMATION 

In Newtonian theory, the so-called Cowling approximation corresponds to neglecting perturbations in the gravitational po- 
tential. The accuracy of the approximation depends on the problem under consideration. For example, while it is not a very 
good approximation for the fundamental /-mode oscillations of a fluid star it can be excellent for the class of gravity g-modes 
which are located in the surface region. Basically, the outcome depends on whether the fluid dynamics induces significant 
variations in the gravitational field. Another problem where the Cowling approximation is useful is that of inertial modes of 
rotating stars. In this case the dynamics is dominated by the Coriolis force and variations in the gravitational potential enters 
at higher orders in the slow-rotation expansion. 

In this paper we make use of the analogous approximation in the framework of general relativity. This problem is 
obviously somewhat more complex owing to the fact that in general relativity the gravitational field is a dynamical entity 
describe d by a tensor field, not just a scalar potential. The relativistic Cowling approximation has a long history dating 
back to iMcDermott et all lll983h who considered /, p and q-mode s in neutron stars. Their appro ach was to simply take 
the equations that had been derived by Thorne and co-workers, e.g. iThorne fc Campolattarol l)l96^ . and neglect the metric 
perturbations. Since the original equations are stated in a specific gauge it is no t clear that this is equivalent to "throwing 
away the gravitational degrees of freedom" . This issue was duly noted by iFinnl l|l98£h who reconsidered the problem and 
suggested a differ ent form of the approximatio n. However, even then the gauge problem was not considered in detail. It was 
later suggested bv lLindblom fc Splinterl il99fj) that the original formulation gave better results (compared to fully numerical 
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work) than the modified version for the low order p-modes. One has to treat this comparison with some caution, however, 
because Finn's argument is mainly relevant for the p-modes. 

In this paper we consider torsional oscillations in the crust (as opposed to the polar perturbation case discussed above). 
We shall pay special attention to the gauge problem and define our approximation as neglecting the coupling between the 
gravitational and matter degrees of freedom rather than just setting the variations of the gra vitational field to zero. It w ould 
be interesting to consider this approach also in the polar case (e.g. using the formalism of lOerlach fc Seneuptal l)l979h ') in 
order to compare with the work cited above. 



Al Axial perturbations in the decoupling limit 

In this appendix we make contact with recent work on relativistic elasticity theory. By default this means that we will have 
to be somewhat technical. However, as the complete problem formulation and all astrophysically relevant results are provided 
in the main text, the material provided here is only intended for readers that are interested in the technical details. 

The main results (i.e . the p e rturbation equations) that we will describe agree with previous work in the relevant limits 
[see lSchumaker fc Thornel il983l) : lYoshida fc Leel 12002k iMessios et alJ l200ll) : ISotani et alJ foodd) ]. but we extend the results 
in the literature to include, in principle, anisotropic and strained backgrounds. The formalism is also valid, with minor 
alterations, for non-static backgrounds. 

We assume that the full perturbed spacetime is axisymmetric with the Killing vector generating the symmetry denoted by 
rf (we use lowercase Latin letters to represent spacetime indices). Since the background around which we will later linearise 
is assumed to be spherically symmetric this only implies trivial restrictions of generality. We write the full spacetime metric 
as 

g a b =-Ub +Ffi a fi b , fi a = F~ l ri a , J7 a _L ai) = 0, (Al) 

where F = ffr\ a — r 2 sin 2 9 is the norm of the Killing vector. The tensor field _L a t, is a metric on t he manifold of K illing vector 
flow lines although it is generically only a metric on a submanifold of the background spacetime. IKarlovinil j2002h has shown 
that under these circumstances the axial perturbations with arbitrary matter sources are governed by the gauge invariant set 
of equations 2 

V b (FQ ab ) = K J a , (A2) 
V a J a = 0, (A3) 
where k = 8ttG/c 4 (= 8 7r in geometric units) is the coupling constant in Einstein's equations, 

Q ab = 2V [a 8n b] , (A4) 

encodes the metric perturbations and 

J a = 25(± ab n c T bc ), (A5) 

where T bc is the stress-energy tensor, describes m atter perturbatio ns. In this section the symbol S is used to denote a perturbed 
quantity in an arbitrary gauge. As discussed bv IKarlovinil (|2002l ) the axial gauge transformations are generated by a vector 
field £ a given by 

C = fav a , r ] a Vafa = 0, (A6) 

leading to 

8/la -> Sfl a +V a f G . (A7) 

If we assume that the metric degrees of freedom are weakly coupled to matter (which should be a reasonable approximation 
in the relatively tenuous crust of a neutron star) we may make the approximation that k — + in the perturbation equations 
(we obviously still keep the gravitational constant for the background). We are then left with decoupled equations for the 
gravitational waves and the matter current J a . 

Note that the gauge invariant two-form Q ab is given solely by the geometric perturbations. Hence, to the extent that the 
coupling may be ignored, the gravitational-wave degrees of freedom (e.g. the w-modes of the star) can be described in a gauge 
in variant manner. We refe r to equation IA2I with the right hand side set to zero as the "inverse Cowling" approximation, 
cf. lAndersson et alJ lll996l) . In order to clarify the content of this equation, let us write it out in Regge- Wheeler coordinates. 
Assuming a harmonic time dependence and separating out the angular dependence in the standard manner we arrive at 



2 



2 



Z = VZ, (A8) 



2 The result is, in fact, more general than this, see IKarlovinil J2002fl for details. 
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where 

_ _ n l fdr\ 2 ld 2 r 2v (l- l)(/ + 2) e 2v r k, , „ , .1 

^VUw "rd^| + e H " = — [Ki + l)r+2(p-Pt)-6m(r)J, (A9) 

and 

4^ = e" - \ (A10) 
dr* 

For matter that is unable to support strain (s uch as perfect fluids) this is just the standard result for axial ui-modes 
Jchandrasekhar fc Ferrari Il99ll lKokkotaslll994) . If, on the other hand, the matter is able to support stresses, Eq. 1A8I 
generalises the ui-mode problem in the inverse Cowling approximation to arbitrary matter. Note that the shear modulus 
enters only via the necessary distinction of pressures, and hence via the background equations. Any direct appearance is 
removed by the inverse Cowling approximation. 

The matter current J a on the other hand depends on the perturbed metric and is therefore not as nicely decoupled. 
Having defined the "gauge invariant" 3 inverse Cowling approximation it is natural to demand that Q ab = in the Cowling 
case. Because of HA7t it is clear that this amounts to letting 8fi a be at most a gradient of a scalar (gauge) function. The matter 
equations will clearly depend on the type of matter. Here we will only consider conformally deforming perfectly elastic matter 
jKarlovini fc Samuelssonlboosl : Icarter fc Quint,analll972F) . although the generalisation to arbitrary matter is straightforward. 



A2 Elastic matter 

For elastic matter the current takes the form l)Karlovini fc Samuelsson1l200fih 

J a = 2(p + P t)FS ab V b (5$-f G ), (All) 

where pt is the pressure in the directions tangential to the surfaces of sphe rical symmetry, &4> is the perturbed (^-coordinate on a 
referen ce space keeping track of the relaxed matter configuration [see e.g. iKarlovini fc SamuelssorJ ( 200.^ ; ICarter fc Quintanal 



( 1974)] and S ah is a metric tensor field that describes axial shear-wave propagation, 

S%=diag(-l,t;?,i;?,0), (A12) 

where v 2 and v 2 are the speed of shear waves in the radial and tangential direction, polarised in the orthogonal tangential 
direction, respectively. In the isotropic limit these speeds are equal and are given by 

v 2 = (A13) 

P + P 

wher e fi is the shear modulus. Under a gauge transformation generated by the vector field (IA6II we have (IKarlovini fc SamuelssorJ 
2006) 



54>^64> + f G , (A14) 

which means that, in any gauge, we have 

J a = 20 + Pt )FS ab V b 84>, (A15) 

whenever Q ab = 0. Therefore (IA3t is also "gauge invariant". 

An intuitive approach to the relativistic Cowling approximation is to start out with the conservation equation for the 
stress-energy tensor and drop the perturbed metric, obtaining 

5{V a T ab ) = V a ST ab = 0. (A16) 

From the results of lKarlovini &: Samuelssonl ll2006l) it is easy to see that the axial part of the perturbed stress-energy tensor, 
dropping the perturbed metric, is given by 

ST ab = 2(p + pt)r) (a S b)c V c 5$. (A17) 

It is now straightforward to convince oneself that the equation resulting from 1A3I is identical to the projection of IA16t 
along rf . All other components vanish identically. Hence, our geometrical way of deriving the Cowling approximation agrees 
with the intuitive approach. The advantage we have gained is that we now have a clear picture of how gauge issues arise. 
Most importantly, we have qualified the approximation as neglecting the coupling between geometrical degrees of freedom and 
matter degrees of freedom. This strategy could prove useful in other contexts, e.g. in a derivation of the Cowling approximation 
for polar perturbations. 

3 The quotation marks should be taken as a reminder that, since we are not using the full linearised equations, the solutions that we 
obtain do not consistently approximate a one-parameter solution to the full Einstein equations to the prescribed order. It is therefore 
somewhat dubious to discuss gauge invariant identifications of points in the perturbed spacetime to the corresponding points in the 
background. 
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We now proceed to derive explicit equations for our elastic crust. We start by separating out the angular dependence. 
Assuming also a har monic time depen dence we write 8<j> — e C(9)F(r), where C(9) turns out to be given by a G egen- 

T in 



bauer polynom i al, see iKarlovii rl <|2002 |). The function F has the s ame meaning a s Y in ISchumaker fc Thornel il98. 
lYoshida fc Leel l)2002h . Wlr in iMcDermott et, alJ l|l988h and F of iMessios et alJ (|200lh if the relevant limits are taken. In 
Schwarzschild coordinates, denoting derivatives with respect to the Schwarzschild radius r by a prime, we find 

F" + A'F' + BF = 0, (A18) 

where 

e A = r 4 e v -\p + Pt )v 2 r , (A19) 

«?(I-DC! + 2)1 (A2Q) 



B= —, 

These equations need to be complemented by boundary conditions at the top and bottom of the crust region. These 
conditions are naturally taken to be that the traction t a = r b 8T a b vanishes at the boundaries. We then have 

t a = (p + pt)r) a vlr b V h 5j> = 0. (A21) 

That the traction should vanish follows from the fact that the shear modulus, and hence the shear wave speed, is zero in the 
fluid or vacuum surrounding the solid. In terms of F we get 



t = r 2 F' = 0, 

at the boundaries. We now define a rescaled traction according to T = e / 
order system, 

T' = 

F' = e~ A j 



-e A BF 



(A22) 

2 t = e A F' and use this to reduce l|A18(l to a first 

(A23) 
(A24) 

This system is convenient for numerical integration since there are no derivatives of the equation of state parameters such as 
the shear modulus and the density (which are only known in tabular form). When the background is isotropic, i.e., when the 
two shear-wave speeds are identical v r = v t = v, the system simplifies considerably and one obtains 

(A25) 

B = < " ■ — - — — . (A26) 



-A - 



(J-l)(J + 2) 



The boundary conditions are now simply that T = at the top and bottom of the crust. 

One can check that this s ystem is equivalent to that obtained bvlYoshida fc Leel ll2002l). which has t he correct Newtonian 
limit JlVIcDermott et al1ll988h . The equations of lSchumaker fc Thornel dl983h and IMessios et alJ J200lt) also reduce to these 
equations when the metric perturbations or magnetic fields are set to zero, respectively. 



APPENDIX B: THE CRUST THICKNESS 

In order to express the estimated crust mode frequencies in terms of just M and R we need to find the crust thickness as a 
function of these variables. To this end we approximate both the mass function and the metric function A as being constant 
and given by e" 2A = 1-2/3, where f3 = M/R. We also use the fact that the pressure is negligible compared to p and in the 
crust and take the equation of state to be a simple polytrope, 

P = kp^ T , (Bl) 

with k and V constant. With these approximations the equation of hydrostatic equilibrium becomes 

' M 2X /-don 

p » -p-je , (B2) 

which may be solved to yield 

XP 1/x «fce a ^+C, (B3) 

where C is an integration constant and x = r/(T — 1). Denoting the pressure at the crust-core interface (at r = R c ) by p c 
(and defining a corresponding transition density p c ) we solve for C to find 

C = a~^-e 2 \ (B4) 

tic 

where a = XPc/Pc- Putting this into HB3t and setting p(R) = at the surface we obtain 



e2XM iR-R-J +a ~°> W 



1 1 



Neutron Star Astero seismology 13 



Numerical data 




P 

Figure Bl. Fit of A/R with a = 0.02326 to numerical results for a large sample of stellar models. Note especially that A/R is (very 
nearly) a function of fi only, and that the simple one-parameter fit suggested by the polytropic approximation is quite accurate. 



A - ^e 2X + l) 1 . (B6) 



solving for A/R — (R — R c )/R we finally obtain 

R " \a 

We see that this is a function of the compactness /3 only, something we have confirmed in the numerical study, see Figure iBll 
The value of a obviously depends on the equation of state. For our tabulated equation of state it can be written 

a^ 0.0047x eff , (B7) 

where Xcff is an effective value. For aT = 4/3 polytrope we obtain a ~ 0.019. We also fitted the expression HB6t using a as a 
free parameter to the numerically obtained data and found 

a « 0.02326, (B8) 

agreeing well with the estimated value. This value corresponds to x ~ 5 or F w 4. In our detailed numerical fits described in 
the main text we use <B6t . However, if one is satisfied with an accuracy at the 10% level then the simple approximation 

is relevant for all neutron stars (since M/R >2x 10~ 2 ). 
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